Experimental and digital investigations of heterogeneity in lower cretaceous carbonate reservoir using fractal and multifractal concepts

Characterization and prediction of reservoir heterogeneity are crucial for hydrocarbon production. This study applies the multifractal theory using both numerical and experimental data to characterize quantitatively the heterogeneity of pore structures in Lower Cretaceous limestone reservoir from the United Arab Emirates. Fractal dimensions calculated from three dimensional digital images showed good correlation (R2 =  + 0.69) with experimental high-pressure mercury injection (HPMI) measurements. Moreover, both experimental and numerical fractal dimensions correlate well with experimental HPMI porosity measurements. Multifractal parameters such as the non-uniformity degree of the pore structures Δα, the asymmetry degree in the vertical axis Δf(α), the concentration of pore size distribution α0 and the asymmetry degree in the horizontal axis Rd estimated from digital and experimental data correlated well and revealed ability to quantitatively describe samples heterogeneity. The ranges of digital and experimental multifractal parameters provided the means to differentiate between homogeneous and heterogeneous samples.


Geological setting
The Thamama Group (Barriasian-Aptian) is a carbonate reservoir that was deposited in shallow-marine, lowenergy carbonate ramp and contains four formations: Habshan, Lekhwair, Kharaib and Shuaiba.Data from this study comes from Lekhwair Formation specifically from the Upper Thamama Zone D. The main depositional facies of Lekhwair Formation include high-frequency coarsening upward cycles of shallow-subtidal skeletal-Bacinella floatstones, skeletal-peloidal wackestones; mud-dominated packstones, capped by skeletal-ooidal grainstones 37 .The carbonate samples were selected carefully to represent six different reservoir rock textures defined within the Upper Thamama Zone D. The Upper Thamama Zone D frequently is considered as low resistivity pay zone (LRPZ) due to the presence of multimodal pore network 37 .

Samples and experiments
Samples used for this study include six core plugs S 1 -S 6 of 25 mm in diameter from Lower Cretaceous, shallow marine limestone reservoir in the UAE covering different depths and textural characteristics.The selected limestone samples are composed of allochems dominated by skeletal fragments, ooids, peloids, and oncoids (Table 1).
Table 1.Brief description of the different samples, core plug features and their limestone textures.The limestones exhibit multi-modal pore network consisting of macropore (e.g., intergranular, moldic, vuggy pores) and micropores within skeletal fragments and between micrite particles, the distribution of which is controlled by depositional textures and diagenetic alterations 30,37 .X-Ray MCT and NCT are non-invasive acquisition methods used to reconstruct 3D tomographic image models of rock 38 .These systems consist of X-ray emitting source and detector receiving the attenuated X-Ray signal crossing the sample.To generate the 3D image, a series of acquisitions were obtained at different angles.Projected data are reconstructed numerically to generate the 3D grey level image representing the sample.The attenuation of incident X-rays is related to the sample inner density.Thus, high grey levels denote the solid phase whereas low grey levels are pores.To investigate the pore geometries and grain morphologies, a segmentation procedure needs to be implemented to separate solid and porous phases.The samples were scanned at coarse scale (20 μm resolution) using a Zeiss Xradia X-ray microcomputed tomography scanner (Fig. 1).Subsequently, two 5 mm high subsamples representing the same sample texture were extracted physically from each core plug by visual inspection of the 3D images (Fig. 2).
The first subsample was used for HPMI measurements to experimentally characterize the capillary pressure curves.This approach involves introducing mercury into samples at increasing pressure, ranging from 0.5 psi to 60,000 psi.The analysis of mercury intrusion volume at various pressures allows for the derivation of a distribution function for pore throat sizes, along with determining the permeability and porosity of the subsample.This experimental method is valuable for rock typing, analysis and interpretation of core data.The results of the pore-throat distributions have yielded suitable information, allowing for the selection of the optimal scanning resolution that effectively captures the majority of pore features.The best resolution that X-ray MCT systems can reach depends on the energy and detector size.As the field of view is constrained by the detector size, smaller subsets were extracted physically from subsamples and scanned at higher resolutions to capture most pore features based on HPMI experimental results (Fig. 3).Furthermore, when experimental results reveal that most pore throats were below the micron, then NCT was implemented to capture pore network.For instance, the 0.5 μm resolution scan of sample S 5 revealed few inter-particle pores that were not connected at this resolution (Fig. 4).Furthermore, the pore throat distribution curve of the same sample confirmed the presence of most pores below the 0.5 μm resolution (Fig. 5).Therefore, smaller representative subsets were further extracted and  scanned using NCT at a resolution of 60 ηm per voxel as illustrated in Fig. 4h.Table 2 summarizes dimensions and resolutions of X-Ray MCT and NCT scanned data for the six studied samples.

Application of fractal and multifractal theories for HPMI data
Application of fractal theory for HPMI HPMI measurements are extremely valuable because they: (i) represent valuable experimental information to characterize the distribution of fluids within the reservoir and rock transport properties, and (ii) provide capillary pressure curves that are commonly used to infer pore size distributions of rock samples.The capillary pressure and saturation relationship depends on several parameters including pore throats size distribution, grain and pore geometry.Depending on the implemented model, the calculated fractal dimensions may reflect specific  where r is the radius to fill a unit of the fractal object, N(r)is the number of the units having a radius r needed to fill the whole fractal object, b is a factor of proportionality and D f is the fractal dimension.The model assumes that a pore is represented by capillary tube with a length of l and a volume equal to πlr 2 .The number of units N(r) can be calculated based on experimental capillary pressure curves measurements.Furthermore, a common assumption implemented in the model considers l independent of r.When the pore structure of the rock is fractal, D f the fractal dimension can be determined by finding the relationship between S Hg the mercury saturation and P c the capillary pressure as the following equation: This fractal dimension reflects the pore size distribution 11 .Moreover, the estimated fractal dimension model revealed high correlations with petrophysical properties of core plugs 11,39 .Other models can provide fractal dimension reflecting the characteristics of pore volume distributions in three dimension.Friesen et al. 20 implemented a model describing measurements of the pore volume of a number of coal and char samples by mercury intrusion porosimetry.The fractal dimension was determined from the relationship in the following equation: where b is a constant.
In addition, some models can provide a fractal dimension characterizing the roughness of the pore surface Zhang and Li 40 .In our study, we focused on the model proposed by Li 21 as one of our goals is to find the correlation between fractal dimension and petrophysical properties such as porosity and permeability in the samples.
Based on Eq. (1) the log − log plot of S Hg and P c was obtained from HPMI experimental measurements for the Sample S 1 as illustrated in Fig. 5a.The pore throat radius distribution revealed in Fig. 5b abimodal behaviour and the cumulative distribution function showed a dramatic increase at pore throat size of r 36% = 10 μm corresponding  to the value log(P c ) = 0.87 PSI at S Hg = 36%.Using (2), the fractal dimension of pores D f was estimated as 2.24.
The plot revealed a poor linear relationship between the variables with a coefficient of determination R 2 = + 0.53.However, the curve showed a double-fractal characteristic depending on the range of the capillary pressure P c . Figure 6a illustrates a strong linear relationship between S Hg and P c with a coefficient of determination R 2 = + 0.98 for all values less than log (P c ) = 0.87 PSI.The same strong relationship is observed for the values larger than log (P c ) = + 0.87 PSI with a coefficient of determination R 2 = + 0.95.Similar result was observed in several previous studies when deriving fractal dimension D f from HPMI experimental measurements Lai and Wang 41 .Moreover, many studies employed multivariate statistical analysis revealed that porosity and permeability showed better results when correlated to pore radius r p where p is the percentage of saturation of the non-wetting phase 42,43 .Therefore, the Swanson's method was implemented to determine the segmentation point for the fractal dimension from each HPMI curve for the studied samples 44 .Several researchers used the position of the maximum value of the plot of S Hg /(P c ) versus S Hg , which is known as the Swanson point, to find the transition between highly and poorly connected pores from HPMI curves 42,43 .The triangular purple point represents the Swanson point of the sample S 1 (Fig. 6b).Based on the Swanson's point, the fractal dimension of the small pores D S and the large pores D L of sample S 1 were estimated from the slopes using Eq. ( 9) respectively as D S = 4.24 and D L = 2.12.This result indicates that the fractal dimension D S = 4.24 of sample S 1 is greater than 3.0 contradicting the Euclidean dimension.Therefore, this range of pores has a low impact on pore structure evaluation.The fractal dimension of large pore throats D L = 2.12 which is conformed to the Euclidean dimension.

Application of multifractal theory for HPMI
The Box-Counting (BC) method is a standard technique used to estimate multifractal parameters revealing the complexity of a data.This method analyses patterns at several length scales to capture self-similarities by zooming in and out in the data.Furthermore, the technique requires the measurements to be equally spaced with a scale of ε.In our study, HPMI measurements were irregularly spaced so linear interpolation were used to obtain regularly spaced data points divided into N = 2 10 = 1024 sub-intervals.Practically, for each scale ε, the mass probability function P k (ε) of the k th interval is defined as the ratio between N i (ε) the pore volume of a k th interval and N t the total pore volume.The mass probability function is an exponential function of scale ε as defined as in the following equation: where α k is the singularity index.
Noted that the number of boxes N(ε) increases exponentially with the scale ε as in the following equation: where N α (ε) represents the number of boxes with singular strength in the interval [α, α + dα] and f(α) is the singularity spectrum 12 .
The generalized dimensions D q and the partition function χ(q,ε) are defined for q > 1 as in the following equation: where χ q, ε = j P q j (ε).For q = 1, D 1 is expressed as in the following equation: The mass exponent τ q and D q the generalized dimension are related through the following equation: For homogeneous objects the mass exponent τ(q) reveals a linear relationship with moments q.Conversely, for heterogeneous objects, the slope of τ(q) may change with respect to q and the deviation is related to the degree of the heterogeneity.
In addition, the multifractal theory provides a relationship between the singularity spectrum f(α) and the singularity α as in as in the following equations: Figure 7 illustrates the mass exponent and singularity spectrum curves for sample S 1 .The variable α 0 represents the concentration of pore size distribution corresponding to the maximum value of the singularity spectrum as illustrated in Fig. 7.The parameter Δα = α min -α max denotes the non-uniformity degree of the pore structures analysed, the larger Δα the higher the data spatial complexity, where α belongs to the interval [α min α max ].The symmetry of the singularity spectrum curve f(α) helps for a quantitative assessment of data heterogeneity.Indeed, uniform data are characterized by symmetric curves whereas heterogeneous ones reveal curve asymmetry 9,45 .The parameter Δf = f(α min ) − f(α max ) denotes the shape characteristics of singularity spectra.The curve f(α) reveals an asymmetry to the right or to the left according to the dominant probability subset.R d = (α 0 -α max )-(α 0 -α min ) represents the asymmetry degree in the horizontal axis of the range [α min , α max ] with respect to α 0 .When the value of R d is positive, it indicates that the porosity distribution is mostly in sparse areas, and the opposite is true for negative values 46 .

Image Segmentation
In order to calculate the fractal and multifractal parameters from 3D X-Ray Digital images, the pore network needs to be identified using an image segmentation method 32 .Indeed, every generated voxel inside the threedimensional X-ray image represents a grey level, coded in 16 bits, associated to the density at this precise spatial location.High and low grey level values denote solid and porous phases, respectively.Standard approaches implement thresholding algorithms to find automatically grey level limits separating these different phases in the image.In this study, we use a common thresholding method called the bi-level segmentation technique.This approach includes the grey levels spatial distribution into the histogram shape information to calculate thresholds separating porous and solid phases.The main advantage of this method is the use of region growing strategy to cope with the fuzzy transition zone representing the unresolved micro porous phase by X-ray imaging.Nevertheless, limitations in image acquisition resolution may lead to an intermediate mode in the histogram related to the presence of micro pores below image resolution.In this paper, the bi-level segmentation method was implemented to find automatically the two thresholds separating the three phases.Lower and higher thresholding grey www.nature.com/scientificreports/level values for the sample S 1 were estimated respectively as 28,253 and 33,023 (Fig. 8b).In order, to calculate the fractal dimension, the intermediate phase was not considered.Each voxel belonging to this intermediate phase represents mixture of pore and grains below image resolution where geometry cannot be captured.Therefore, the binary image used for fractal dimension estimation implemented the lower thresholding value.

Image fractal and multifractal parameters
The 3D Box-Counting method was implemented to analyse patterns of self-similarities in 3D segmented images by breaking them into a grid of boxes.This approach is used, in general, to approximate the Hausdorff dimension of a fractal dimension 45 .Boxes correspond to cubes used in investigating, by zooming in and out, the geometric complexity in a data at several length scales.In our application, pore structure was studied as main pattern of interest in three-dimensional binary images.The box-counting quantifies the presence of pores in each cube for a specific size of box.Practically, the number of cube boxes is calculated including the pore phase N(ε) for a specific cube side length ε.The procedure was repeated by covering the 3D image by a sequence of boxes of descending cube side lengths ε.Subsequently, the fractal dimension was estimated by the slope of the regression straight line representing the relationship between ln(N(ε)) and ln(1/ε) as in the following equation: The extension of the fractal theory is known as multifractal analysis, and it is used when the geometry of a system cannot be described by a single fractal dimension.Instead, a range of fractal dimensions is defined to fully define these objects.Multifractals are measured using a probability distribution P k in each k th box as P k (ε) ∼ ε α k as in shown in Eq. ( 3), where ε is the box size and α k is the Lipschitz-Holder exponent characterizing the sin- gularity strength in the k th box.Zhang et al. 23 proposed the use of α k factor to measure the level of complexity in the spatial distribution.They suggested that the number of boxes, denoted as N α (ε) , can be used to represent the probability P k of singularity strengths between α and α + dα.Additionally, the size of the boxes, ε, can be related to N α (ε) as following N α (ε) ∼ ε −f (α) where α corresponds to the singularity and f(α) is the singularity spectrum as in (4,9,10).Similarly, the mass exponent τ q , the generalized dimension D q , the partition function χ(q,ε), concentration of pore size distribution α 0 , the non-uniformity Δα = (α min -α max ), the asymmetry degree in the vertical axis Δf = f(α min ) − f(α max ) and the asymmetry degree in the horizontal axis R d of the range [α min , α max ] with respect to α 0 were derived from the same equations in section "Samples and experiments".Figure 8 illustrates the result of image segmentation in addition to the generalized dimensions D q and f(α) the singularity spectrum obtained for sample S 1 .www.nature.com/scientificreports/Numerical simulations of rock properties Simulating permeability at pore scale using Lattice Boltzmann (LB) method on 3D X-ray images involves the following steps: (1) acquisition of 3D images of a porous material, (2) segmentation of the images to extract pore network and solid matrix, (3) assignment of boundary conditions and fluid properties to the pore network, (4) simulation of fluid flow through the pore network, and (5) calculation of permeability from the simulated fluid flow.The method considers the micro-scale features of the porous material, providing a more accurate representation of permeability compared to macro-scale models.In this approach, the fluid is modelled as particles moving in a lattice structure and their movement is described using a time and space distribution function.The governing equation calculates the density f(x,t) of the particles at each iteration, taking into account both the streaming and colliding terms, where x represents the location and t represents the time as in the following equation: The particle movement direction in the lattice is indicated by e i , where i is the index for a specific direction.The particle velocity is u i , τ is the relaxation time controlling the rate of approach to equilibrium, Ω is a collision operator required to satisfy the conservation of total mass and total momentum, and F is an external force term.The LB method has several advantages as it calculates the particle velocity at each iteration using only the velocities of surrounding particles, making it easily implementable on high performance computer clusters.A multi-relaxation-time model based on the Lattice Bhatnagar-Gross-Krook (LBGK) scheme for the 3D lattice was used in this study, following a D3Q19 lattice model 47 .The flow was set to be periodic between the inlet and outlet and the bounce-back rule was used to manage solid-fluid boundaries.Additionally, a no-slip boundary condition was implemented at the solid-fluid interface.When the steady state is reached, the permeability is estimated using the unidirectional Darcy's law as in the following equation: where K is the absolute permeability, ΔP represents the pressure difference along the length of the sample, A is the surface section area, Q stands for the flow rate and µ is the dynamic fluid viscosity.The permeability value estimated is in lattice units and are converted to a real value using the resolution of the scanned image.Moreover, the porosity can be determined by calculating the ratio of the number of voxels within the segmented pore space to the total number of voxels in the image.This value is calculated directly from the segmented 3D images.

Numerical and experimental fractal analysis of pore structures
Comparison between the simulated and experimental properties showed a good agreement for both porosity and permeability rock properties with coefficient of determinations of R 2 = + 0.69 and R 2 = + 0.98, respectively (Fig. 9a,b).The correlation between experimental and simulated permeability values is higher than those of porosity because permeability have larger magnitude variation, from 0.02 to 800 mD.These observations indicate that the three-dimensional images captured successfully digital pore networks representative of the real pore space in subsets used for HMPI experiments.This result confirms that image segmentation and Lattice Boltzmann methods provide accurate results for porosity and permeability when implemented on high quality resolution images as reported in several studies 30,48,49 .Figure 9c reveals a fair agreement between the fractal dimension D L derived from HPMI and the fractal dimension FD L calculated form 3D images with a coefficient of determinations R 2 = + 0.69.This observation suggests the existence of a link between pore structures description using image fractals and experimental rock properties in the studied carbonate samples.
Furthermore, the image fractal dimension FD L showed a positive correlation with porosity with a coefficient of determination R 2 = 0.56 for a direct linear fitting relationship (Fig. 9d).This observation suggests that high porosities correspond to high image fractal dimension FD L values because high porosity images usually have more compact pore structure.Figure 9e revealed no correlation (R 2 = + 0.13) between the image fractal dimension and the experimental permeability measurements obtained from HPMI.Indeed, the fractal dimension FD L also called capacity dimension usually captures the average presence of pores rather than their connectivity 29 .Furthermore, a positive correlation was revealed between the fractal dimension FD L and the radius at Swanson's segmentation point with a relatively fair coefficient of determination R 2 = + 0.45 (Fig. 9f).

Heterogeneity characterization using HPMI experiments
Various types of pores with varying sizes develop in carbonate reservoirs, leading to heterogeneity in the structure of the pores (Chen et al., 2016).A brief description of the six samples features and their limestone textures is provided in Table 1.The permeability of samples S 1 and S 2 derived from HPMI results analysis are 401.63mD and 112.84 mD respectively (Table 3).These high values of permeability shows that samples are dominated by large pore throats.This observation is in agreement with the wide bimodal pore throat size distributions with a major peak at around 16 µm (Fig. 3).A second peak appears around 0.113 µm and 0.21 µm for samples S 1 and S 2 , respectively.Both samples show well connected large pores in their macro-porous phases and 3D-MCT scan data, and the presence of inter-particle and intra-particle pores.In contrast, the samples S 3 to S 6 showed low permeability values spanning from 0.03 mD to 6.75 mD, suggesting that the samples are dominated by small pore throats (Table 3).Samples S 3 and S 4 showed bimodal pore throat size distributions, narrower than samples S 1 and S 2 , with a major peak at around 1.2 µm (Fig. 3).A smaller peak appeared for both samples at 0.58 µm as well.
The 3D-MCT image of sample S 3 acquired at 0.5 µm resolution showed a well-connected macropore system and presence of dolomite rhombs within the micrite matrix.The 3D-MCT scan of sample S 4 at 0.5 µm resolution showed that the micro-porous phase exits inside the leached grains in grain supported texture.The HPMI data of samples S 5 and S 6 revealed unimodal pore throat size distributions with a major peak at 0.27 µm and 0.14 µm, respectively.Both samples display relatively narrow distributions by comparison to the other four samples (Fig. 3).The NCT image acquired at 60 ηm resolution revealed that most of the microporous phase exists inside the mud matrix for sample S 5 and exists within the micrite between the oncoidial allochems for sample S 6 .Furthermore, the fractal characteristics of the pore space for the six samples were investigated using experimental HPMI measurements.In percolation theory suggests that the flow properties are predominantly influenced by a characteristic length, which plays a pivotal role in governing the fluid flow in reservoir rocks 42 .An increased fractal dimension indicates a shift in the regular pore morphology towards a more complex form 15,21,41 .Consequently, this transformation results in a decrease in both fluid flow and permeability.Several researchers have used fractal dimension on tight sandstone to predict permeability.For example, Lai and Wang 41 , applied Li's 21 model to assess the fractal characteristics of tight sandstones, revealing the presence of two distinct fractal regions in all tight sandstone samples in their investigation.The fractal dimension value acts as a quantitative measure of reservoir heterogeneity, demonstrating that an elevation in the fractal dimension corresponds to a more complex pore structure, ultimately resulting in decreased permeability.Nevertheless, in our study Fig.  www.nature.com/scientificreports/reveals the existence of a relatively strong direct linear relationship between the fractal dimension and the logarithm of the simulated permeability with R 2 = 0.72.This outcome is consistent with the result reported by Xin et al. 26 revealing an increase of permeability with respect to the increase of fractal dimension in a fractured-vuggy carbonate reservoir.However, Zhang et al. 15 studied fractal dimension of pore structures of Lower Carboniferous carbonate reservoir and showed the existence of an inverse linear relationship with permeability (R 2 = 0.75).
Figure 11 shows the log − log plot of S Hg and P c samples S 2 to S 6 .As for sample S 1 the plot showed poor linear relationships fitting between the variables.Nevertheless, Fig. 12a revealed that all curves have a double-fractal behaviour characterized by the Swanson segmentation points.The curves obtained by plotting S Hg versus S Hg /P c denoted sharp apex for each sample (Fig. 12b).Table 4 summarizes the maxima values obtained for each sample at the corresponding saturation percentage S Hg .The double-fractal characteristic is revealed by piecewise linear regressions showing good fits with a positive coefficient of determination R 2 spanning in the interval [0.87, 0.95] for the six samples.The fractal dimensions of the large pores, denoted D L , ranges from 2.04 to 2.23.In addition, a good fit is revealed for the left segments corresponding to smaller pores in Fig. 12a, the estimated fractal dimensions D S values are all greater than 3.0, which contradicts the Euclidean dimension revealing absence of fractal characteristics.Among, our studied six samples, only sample S 2 may reveal three sections for the log(Pc) versus log(Shg) curves.However, the estimation of the fractal dimension in the intermediary section provided also a fractal dimension larger than 3.0 contradicting the Euclidean dimension.Therefore, we used only D L revealing the fractal behaviour for macro-pores.Moreover, several multifractal parameters were calculated based on the HPMI measurements of the six samples.The mass exponent τ(q) and the singularity spectra f(α) were calculated for the six samples.Figure 13a illustrates slope variations of the mass exponent τ with respect to q for the six samples.Samples S 1 and S 2 show higher changes in slopes values compared to the other four samples.These large variations are due to the higher heterogeneity in pore throat size distributions in these two samples described in both HPMI data and observed in 3D X-Ray images.This observation confirms the ability of multifractal parameters to assess quantitatively heterogeneity in HPMI measurements.The singularity curves for the six samples revealed convex parabolic shape as illustrated in Fig. 13b.The samples S 1 , S 4 and S 6 revealed singularity spectra curves f(α) with a left asymmetry around α equals to 1. Singularity curves of samples S 1 , S 4 and S 6 showed wider left portions with sharper slopes than the right ones.However, asymmetric singularity curves of samples S 2 , S 3 and S 5 revealed wider right

Multifractal characteristics of pore structures using image analysis
The analysis of 3D-MCT images of samples S The singularity curves f(α) for the six samples showed right sided asymmetric convex parabolic shapes as illustrated in Fig. 15.Also, Table 6 summarizes the image multifractal parameters results for the six samples.The concentration of pore size distribution values α 0 ranged between 2.977 and 2.988.Larger α 0 values indicate smaller concentrated distributions.Samples S 1 and S 2 revealed largest Δα values among the six samples respectively 0.413 and 0.416 whereas Δα was in the range [0.357, 0.392] for the other four samples.This observation confirms that the non-uniformity degree of the pore structures Δα captures heterogeneity from 3-dimensional image data.Moreover, Δf(α) values were negative for the six samples indicating that a small probability subset was dominant in these samples.

Correlation between HPMI and digital multifractal parameters
The correlations between multifractal parameters derived from HPMI experimental measurements and three dimensional images were analysed.Figure 16a revealed the existence of a strong linear relationship between the non-uniformity degree Δα calculated from HPMI and images with a relatively high determination coefficient R 2 = + 0.73.Furthermore, the concentration of pore size distribution α 0 showed also a relatively good linear relationship with a determination coefficient R 2 = + 0.69 (Fig. 16b).However, the relationship between the shape characteristics of singularity spectrums Δf(α) derived from images and HPMI experiments showed a weaker correlation with a determination coefficient R 2 = + 0.51 (Fig. 16c).These results show that the multifractality notion could describe independently both digital and experimental representations of the same data and that a correlation exists between these two representations.
The results obtained in this study indicate the relevance of multifractal parameters to describe heterogeneity quantitatively in pore structures from both experimental and digital data.The comparison between multifractal parameters ranges allowed discriminating homogeneous from heterogeneous samples.Nevertheless, this method suffers of lack of clear established reference ranges.In future studies, it is imperative to explore a broader spectrum of experimental and digital measurements to establish statistically more comprehensive parameter ranges for Δα, α 0 , R d , and Δf(α) in carbonate rocks.These parameters serve as quantitative indicators for classifying the degree of heterogeneity in carbonate formations.Subsequently, in a secondary phase of the research, the goal is to identify the most significant multifractal parameters among these, influencing key rock properties such as porosity, permeability, and elasticity.Ultimately, our aim will be to integrate the most influential parameters among Δα, Δf(α), R d , and α 0 as quantitative heterogeneity indicators into rock physics models for carbonate formations.

Figure 1 .
Figure 1.HPMI experiments subsample selection and image acquisition at high resolution for Sample S 1 : (a) core plug photograph, (b) X-Ray 3D-MCT at 20 μm resolution and selection by visual inspection of subsamples, (c) Subsample (trim) to be examined at high resolution, and (d) Subsample used for HPMI measurements.

Figure 2 .
Figure 2. Subset selection and image acquisition at high resolution for Sample S 1 : (a) Location of 4 mm 3 extracted subset on 3D X-Ray MCT at 20 μm resolution, (b) subset location on real core plug, (c) Subset extraction, and (d) Subset scanned at 4 μm resolution.

Figure 3 .
Figure 3. Semi logarithmic plot of experimental HPMI Pore throat distributions of the six samples.

( 1 )Figure 5 .
Figure 5. (a) Log-log plot relating S Hg the mercury saturation to the capillary pressure P c for the sample S 1 , and (b) Plot of pore throat radius distribution and its cumulative function.

Figure 6 .
Figure 6.(a) Log-log plot relating S Hg the mercury saturation to the capillary pressure P c for the sample S 1 illustrates the double-fractal characteristics, and (b) Plot of log S Hg P c versus S Hg illustrating the position of Swanson's point: pore throat radius r 36% = 10 μm for the sample S 1 .

Figure 8 .
Figure 8. Fractal dimension calculation process based on 3D X-ray micro-computed tomography for sample S 1 .(a) The original 3D image, (b) histogram of grey levels with the thresholds implemented for 3D image segmentation, (c) segmented image, (d) the generalized dimension D q , (e) Fractal dimension D = 2.62 calculated as slope of the line ln(N(ε)) as a function of ln(1/ε), and (f) Singularity spectrum f(α) calculated from 3D segmented image.

Figure 9 .
Figure 9. Correlation between experimental and simulated rock properties for the 6 samples: (a) porosity, (b) permeability, (c) experimental and simulated fractal dimensions D L and FD L .Correlations between image fractal dimensions FD L and (d) porosity, (e) permeability, and (f) Radius at Swanson's segmentation point.

Figure 10 .
Figure 10.Relationship between fractal dimensions FD L and simulated permeability log(K) for sample S 1 to S 6 .

Figure 11 .
Figure 11.Log-log plot relating S Hg the mercury saturation to the capillary pressure P c for sample S 2 to S 6 .

Figure 12 .
Figure 12.(a) Swanson's segmentation points for samples S 1 to S 6 , and (b) Log − log plot of S Hg and P c obtained from HPMI experimental measurements for samples S 1 to S 6 revealing double-fractal behaviour.

Figure 14 .
Figure 14.Fractal dimensions calculated as slope of the line ln(N(ε)) as a function of ln(1/ε) in sample S 1 : in red D = 2.62 including all pores and in blue D = 2.27 including all pores with r > 10 µm.

Table 2 .
Subset dimensions and corresponding 3D-MCT and NCT image acquisition resolutions.

Table 3 .
Experimental and numerical porosity and permeability values for the six samples.

Table 4 .
HPMI fractal dimensions, image fractal dimensions and Swanson's segmentation points for the six samples where x represents the saturation of Hg in %.

Table 5 .
Multifractal parameters from HPMI measurements.alargeprobability subset was dominant in these samples.The samples S 2 , S 3 and S5show Δf values that are negative indicating that a small probability subset was dominant in these samples.
Vol:.(1234567890) Scientific Reports | (2023) 13:20306 | https://doi.org/10.1038/s41598-023-47681-wwww.nature.com/scientificreports/ that 23-S 4 produced fractal dimensions values D ranging from 2.53 to 2.64, whereas NCT images gave fractal dimensions 2.59 and 2.74 for samples S 5 and S 6 , respectively.The image fractal dimensions D values (Table4) reveal larger values when compared to the range of experimental fractal dimension obtained from HPMI which are in the range[2.10,2.23].This result was expected as fractal dimension calculation was based on all pores included in images.Therefore, to have a more representative comparison, image fractal dimensions (FD L ) were recalculated keeping only pores with sizes larger than the Swanson's segmentation points.Figure14illustrates the comparison between fractal dimension values calculated including all pores in red and including only pores larger than 10 µm in blue.The recalculated fractal dimension values vary in the range [2.05, 2.38] (Table4).The fractal scaling law, as observed in various studies, suggests that an increase in porous area (in 2D images) or volume (in 3D images) will result in higher fractal dimensions23.Consequently, the decrease of fractal dimension values is in agreement with results reported in these previous studies.Furthermore, samples S 1 and S 2 revealed image fractal dimensions FD L values of respectively 2.27 and 2.38 larger than the four other samples with values ranging in the interval [2.05, 2.28].This observation indicates the effectiveness of fractal dimensions in capturing complexity in pore structures of samples through the imaged pore network.